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1.  INTRODUCTION 


The  Coordinate  Systems  Conversion  Software  (CSCS)  is  a  collection  of  general  purpose  subroutines  and 
functions  which  may  be  used  to  transform  vector  quantities  from  one  coordinate  system  to  another.  It 
is  intended  for  use  by  persons  with  a  good  working  knowledge  of  the  coordinate  systems  and  parameters 
of  interest.  The  routines  are  low  level  in  the  sense  that  a  single  call  is  only  a  single  step  in  the 
conversion  process.  By  combining  several  steps,  however,  the  user  will  find  the  CSCS  to  be  a  powerful, 
general  purpose  tool  which  is  capable  of  solving  a  wide  variety  of  problems. 

The  CSCS  collection  contains  routines  which  fall  into  five  categories:  quaternion  conversion,  rotation 
matrix  construction,  attitude  manipulation,  M50  to  Aries  true-of-date  conversion,  and  vector 
manipulation.  Some  routines  may  foil  entirely  within  a  single  category.  Others  span  multiple  categories, 
thus  forming  a  bridge  between  them.  By  using  such  bridges  the  user  may,  for  example,  convert  a 
quaternion  into  a  rotation  matrix,  which  is  then  used  to  manipulate  a  vector.  The  following  sections 
describe  the  routines  in  the  collection  and  give  several  example  programs  to  demonstrate  their 
application. 

Attitude  sequences,  quaternions,  and  rotation  matrices  are  used  by  the  space  community  to  store  and 
transmit  coordinate  system  conversion  information.  Each  method  has  its  own  distinct  advantages  and 
disadvantages,  and  the  adoption  of  one  method  in  preference  to  another  takes  such  considerations  into 
account.  In  theory,  attitude  sequences  contain  all  of  the  conversion  information  in  the  form  of  three 
rotation  angles.  They  are  easy  to  understand  in  human  terms,  and  translate  easily  into  motions  such  as 
turning  one’s  body,  looking  up,  and  tilting  one’s  head.  If  the  sequence  of  rotations  is  not  known, 
however,  their  meaning  is  lost.  In  addition,  they  must  be  converted  using  trigonometric  functions  before 
they  can  be  used  for  vector  analysis.  Quaternions  store  the  conversion  information  as  four  unambiguous 
numbers.  They  may  be  converted  to  a  form  suitable  for  vector  analysis  somewhat  more  easily  than 
attitude  sequences.  Thus  they  are  favored  in  situations  where  speed  and  machine  storage  considerations 
dominate.  They  are  not  easy  to  understand  in  human  terms  and  cannot  be  used  directly  in  vector 
computations.  Rotation  matrices  consist  of  a  3x3  array  of  nine  numbers.  They  can  be  used  directly  in 
vector  multiplication  and  thus  are  preferred  for  data  analysis  using  computing  machines.  The  CSCS  is 
capable  of  working  with  data  in  any  of  these  forms  and  converting  it  to  any  of  the  other  forms  to  meet 
the  needs  of  the  user. 

In  selecting  a  reference  ECI  coordinate  system  the  space  community  has  commonly  employed  two 
different  systems,  the  MSO  and  the  Aries  true-of-date.  A  particular  mission  will  generally  maintain 
orbital  and  ephemeris  information  in  one  reference  system  but  not  the  other.  The  CSCS  provides  the 
ability  to  convert  data  to  the  other  ECI  system  for  any  date  and  time  selected  by  the  user.  The  M2000 
system  is  likely  to  become  more  prevalent  in  the  future,  and  the  MSO  conversion  procedures  may  easily 
be  modified  for  application  to  M2000. 

Finally,  the  CSCS  contains  vector  manipulation  routines.  One  routine  allows  the  user  to  convert  vectors 
from  one  coordinate  system  to  another.  A  different  routine  gives  the  user  the  ability  to  reverse  the 
process;  that  is,  to  align  an  axis  with  a  vector  and  obtain  the  associated  rotation  matrix.  With  this 
capability  the  CSCS  becomes  a  powerful  tool  which  lets  the  user  construct  coordinate  systems  of  interest 
together  with  their  related  inter-conversion. 


1 


The  Coordinate  System  Conversion  Software  described  in  this  document  was  developed  for  IBM  PC 
compatible  systems  using  Microsoft  FORTRAN.  All  floating  point  variables  and  calculations  are  double 
precision  to  ensure  high  accuracy.  Constants,  such  as  pi,  are  defined  by  the  operating  environment 
whenever  possible.  ANSI  FORTRAN-77  standards  were  adhered  to  so  that  the  code  may  be  transported 
to  other  environments  with  a  minimum  of  conversion  effort. 


2.  NOTES  ON  THE  USE  OF  QUATERNIONS 


Two  references  were  used  in  the  development  of  the  software  contained  in  this  document.  The  NASA 
PATH  documentation  [Cooper,  et  al.,  1985]  was  referred  to  for  all  routines  which  use  quaternions. 
Escobal  [1965]  was  referred  to  for  routines  which  involve  attitude  conversions  and  coordinate  system 
rotations.  The  Appendices  of  Escobal  were  especially  helpful  in  this  respect.  Escobal  also  discusses 
quaternions  in  great  detail  and  relates  them  to  these  other  areas.  NASA  and  Escobal  have  defined  the 
order  of  operations  differently,  however,  and  the  unwary  user  will  quickly  foil  victim  to  this  unfortunate 
circumstance.  The  NASA  definition,  V  ■  q  V  q*,  works  just  as  well  as  that  of  Escobal,  V  -  q*  V  q, 
as  long  as  one  does  not  mix  conventions. 

In  the  development  of  the  quaternion  routines  presented  here  the  NASA  convention  was  followed. 
Much  of  the  material  discussed  by  Escobal  is  still  relevant,  but  the  derivations  look  quite  different  in  the 
end.  If  an  Escobal  convention  quaternion  is  encountered,  it  is  recommended  that  the  user  convert  it  to 
its  conjugate  before  proceeding.  This  may  be  done  with  the  QTPOSE  routine,  which  will  return  a 
quaternion  consistent  with  the  NASA  convention. 


3.  DERIVATIONS  OF  COORDINATE  SYSTEM  TRANSFORMATIONS 


We  introduce  the  quaternion,  f,  in  terms  of  four  parameters  which  are  components  of  q  as  follows: 

f  -  (fof1.f2.J3)*  (1) 


We  also  define  the  multiplication  of  two  quaternions,  r  =  (y  [HoUinger,  1989] 


(2) 


It  is  sometimes  more  convenient  to  think  of  quaternions  as  being  comprised  of  scalar  and  vector 
components.  With  scalar  component  qQ  and  vector  components  q{\+q$+q$K  we  define  the  quaternion 
f  [Escobal,  1965]  as 


2 


where  i,  j,  and  k  are  the  hyperimaginary  numbers  satisfying  the  conditions 


i2  -  j2  m  k1  «  -1 

a  -  -j*  -  *  (4) 

jk  -  -kj  m  i 
ki  ■  -tk  -  j. 

Using  this  definition,  it  is  now  possible  to  carry  out  quaternion  multiplication  in  the  same  manner  as 
polynomial  multiplication.  We  must  take  care,  howevet,  to  preserve  the  order  of  operations,  since  the 
hyperimaginary  numbers  clearly  do  not  commute. 

We  now  define  the  conjugate  of  a  quaternion  as 

f*  ■  (5) 


We  note  that  a  vector  in  three-dimensional  space  with  components  qxi+  q3k  may  be  expressed  as 
a  quaternion  with  a  scalar  component  of  zero  in  quaternion  space 

f  *  O  +  fyl  +  fj/  +  ?3*.  (<0 


3.1  CONVERSION  OF  A  QUATERNION  INTO  A  ROTATION  MATRIX 

The  transformation  of  a  vector  Vto  V  is  accomplished  by  the  quaternion  multiplication  [Cooper,  et 
ol.,  1985] 

Y  »  f  Vi*.  W 

This  is  equivalent  to  multiplication  of  Vby  the  rotation  matrix  R 

Y  »  RV.  W 

Carrying  out  the  quaternion  multiplication, 

y  -  (?o  +  ?i*'+  V  +  ?3*)  (v  +  V  +  v3*)  (9) 

*  we  reduce  this  expression  to 
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(10) 


*  [2 (ff0?3  +  ?l?2)vl  +  (*o“  f ?  +  h  “  *3 )  W2  +  2(*2*3  -  9Q?!  )  V3  ]/ 

♦  [2(*1*3  -  *)*2  )  W1  +  2(*0*1  +  *2*3  )  v2  *  (?0  ‘  *1  '  *2  +  *3  )  tf3  ]*• 


The  elements  of  J?  are  easily  identified  as 
[2222 

ffo  +  *i-*2-*3 


R 


2(90ff3+«l«2) 

2(9i93-90?2) 


2(9i*2  “*0*3) 
2222 

«0-?l+?2-«3 

2(f0?l  +  *2*3) 


2(96*2  +  lift) 

2(92*3“  *0*1  )  • 

2  2  2  2 
*0-*l-*2  +  »3 


(ID 


The  transformation  from  the  prime  to  the  unprimed  system  is  accomplished  using  the  transpose  of  R, 

V-Rr\r.  (12> 


The  quaternion  equivalent  is 

V-9  *Vq.  <13) 

We  see  that  talcing  the  quaternion  conjugate  is  the  equivalent  of  matrix  transposition,  and  use  these  terms 
interchangeably. 


3.2  CONVERSION  OF  A  ROTATION  MATRIX  TO  A  QUATERNION 

Using  the  quaternion  expressions  for  the  elements  of  our  rotation  matrix  derived  above,  we  establish 

^11  +  ^22  +  ^33  + 1  M  4*o 

^11 -y?22-J?33  +  1  *  4*f  (14) 

“  ^11  +  ^22  “  ^33  +  1  "  4  *2 
■^11  ~^22  +  ^33  +  1  "  4  *3  * 

Adding  these,  we  show 
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(15) 


2  2  2  2 
*0  +  *l+*2  +  *3 


1, 


» 


a  useful  identity. 


The  relationships  above  give  the  magnitude  of  the  quaternion  components,  but  not  their  signs.  Changing 
the  signs  of  all  four  quaternion  components  does  not  affect  our  rotation  matrix,  so  we  may  make  our 
first  choice  arbitrarily.  If  the  square  root  of  the  expression  for  qQ  above  is  not  zero,  we  select  the  positive 
root.  We  then  solve  for  the  remaining  quaternion  components  by  referring  to  our  rotation  matrix  and 
seeing  that 


Rtt-Rn 


*3 


(16) 


If  qQ  is  zero,  we  use  the  square  root  of  the  expression  for  f,  above  and  proceed  in  a  similar  manner. 


3.3  THE  PRODUCT  OF  TWO  QUATERNIONS 

We  are  interested  in  finding  the  equivalent  quaternion,  f  ',  for  the  transformations  f  followed  by  f '. 
In  terms  of  a  vector,  V,  two  successive  applications  of  the  quaternion  transformation  operator  give 

v”  -  f'(fV)f'*.  (17) 

It  can  be  shown  that  the  hyperimaginary  numbers  are  associative,  thus 

V"  m  (*'t)  V(q* §'*).  (18) 

It  follows  that 
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f '  f  -  (ft '  +  fl'  *‘+  *2  'j+  ft '  *)(ft  +  ft  *+  ft  J  +  ft  *) 
9  (ft  fo  *  ~  ft  ft  ’  ~  ft  ft '  ~  ft  ft  * ) 

♦  («bfr+ft'fiM2'«»"ft'f2)< 

and 

fV*  *  (fo-li*-fci-f3*)(«b'”?r*”f2V~fe'*) 

“  (ftft'-ftft'-ftft'“ftft') 

+  (-ftft'-ft'ft-ft'ft+ft'ft)  * 

+  (-ftft'-ft'ft  +  ft'ft"ft'ft)  / 

+  (-ftft'-ft'ft  +  ftft'-ft'ft)  *  • 

From  these,  we  identify  the  components  of  f': 

ft"  “  ftft'-ftft '-ftft'-ftft' 
ft"  ■  ftft'  +  ftft'-ftft'  +  ftft' 
ft"  “  ftft'  +  ftft'  +  ftft'"ftft' 

*3  “  ft  ft  “  ft  ft  +  ft  ft  +  ft  ft  ’  • 

We  may  now  use  f '  to  operate  on  V  directly 

V"  -  f"  Kf"\ 


(19) 


(20) 


(21) 


(22) 
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3.4  ROTATE  AN  AXIS  TO  BE  COLLINEAR  WITH  A  VECTOR 


This  problem  is  treated  in  a  general  way  by  Escobal  [1965).  A  single  rotation  about  an  arbitrary  rotation 
axis  is  always  sufficient  to  make  one  vector  collinear  with  another.  Defining  a  unit  vector  e  along  the 
rotation  axis,  and  a  rotation  angle  ®, 


R 


2 

COS®  +«j  (1  -cos®) 
ej  e2(l  -  cos® )  -  e3sin® 


«!  e2(l  -  cos®)  +  e3sin®  «j  e3(l  -cos® )  -  e2sin® 
2 

cos®  +  e,(l -cos®)  e2e3(l-cos®)  +  eiSin® 


«je3(l -cos®) +  e2sin®  ej^Cl-cos®)-*^^®  cos®  +  e*  (1  -cos®) 


(23) 


We  simplify  the  problem  by  requiring  that  the  first  vector,  U,  be  equal  in  magnitude  to  the  second  vector, 
V.  Furthermore,  we  require  U  to  point  in  the  direction  of  the  n  axis.  The  angle  ®  is  the  direction 
cosine  between  V  and  the  n  axis: 


®  ■  arccos — - . 

\V\ 


(24) 


The  transformation  V  =  RU  gives  three  equations  with  three  unknowns,  since  t»n  =  /V/  and  the  two 
other  components  of  U  are  each  zero. 


For  example,  if  we  select  U  to  point  along  the  z  axis  (see  illustration),  we  have 


0 

|  V|{e1e3(l-cos®)-e2sin®} 

V  •  R 

0 

m 

1  ^|{«2<3(1  -cos®)  te^m®} 

\V\. 

|  V\  (cos®  +  «3(1  -cos®)} 

Solving  for  e3 


v3 

Im 


-cos® 


/(1-cos®) 


(25) 


(26) 


as  expected,  since  the  axis  of  rotation  must  lie  in  the  xy  plane.  Solutions  for  the  other  components 
become 
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4.  GENERAL  DESCRIPTION 


4.1  QUAT.FOR 


AEQA  This  function  compares  two  attitude  sequences.  A  value  of  .TRUE,  is  returned  if  they 
are  equivalent  and  a  value  of  .FALSE,  if  they  are  not.  Acceptable  error  limits  are 
specified  in  the  logical  function  REQR.  All  angles  are  in  radians.  See  description  of 
subroutine  ATOR  for  details  of  parameters. 

ATOR  This  subroutine  converts  rotation  sequence  about  IAXIS(l),  IAXIS(2),  IAXIS(3)  into 
the  corresponding  nine  element  rotation  matrix.  Allowed  values  are  1- roll* x-axis, 
2-pitch-y-axis,  and  3-yaw-z-axis.  The  same  rotation  axis  cannot  be  used  twice.  All 
attitudes  passed  to  this  routine  must  be  in  radians.  For  example,  IAXIS(l)-2, 
IAXIS(2)-3,  IAXIS(3)«1  specifies  a  pitch,  yaw,  roll  sequence  of  angles  ATT(1), 
ATT(2),  ATT(3)  respectively.  Returns  the  rotation  matrix  R. 

QEQQ  This  function  compares  two  quaternions.  A  value  of  .TRUE,  is  returned  if  they  are 
equivalent  and  a  value  of  .FALSE,  if  they  are  not.  Absolute  errors  on  the  order  of  IE-6 
are  acceptable. 

QMULTQ  This  subroutine  multiplies  two  quaternions,  Q(3)-Q(2)Q(1).  It  is  equivalent  to  two 
successive  rotations,  0(1)  followed  by  Q(2).  Note  that  quaternion  multiplication  is  the 
equivalent  of  matrix  multiplication  R(3)-R(2)R(1)  for  the  same  rotations. 

QTOR  This  subroutine  converts  a  four  element  quaternion  into  a  nine  element  rotation  matrix 

for  coordinate  system  transformations. 

QTPOSE  This  subroutine  transposes  the  quaternion  Q1  into  Q2. 

REQR  This  function  compares  two  rotation  matrices.  A  value  of  .TRUE,  is  returned  if  they  are 

equal  and  a  value  of  .FALSE,  if  they  are  not.  Absolute  errors  on  the  order  of  IE-6  are 
acceptable. 

RMULTR  This  subroutine  multiplies  two  rotation  matrices,  R3  -  (R2)(R1).  It  is  equivalent  to  two 
successive  rotations,  R1  followed  by  R2. 

RMULTV  This  subroutine  multiplies  a  vector  by  a  rotation  matrix,  V2  -  (R)(V1). 


9 


QUAT.FOR  (CONT.) 

RTOA  This  subroutine  converts  a  nine  element  rotation  matrix  into  the  corresponding  rotation 

sequence  about  IAXIS(l),  IAXIS(2),  IAXIS(3)  where  1-roll-x-axis,  2-pitch-y-axis, 
and  3»yaw«z-axis.  The  same  rotation  axis  cannot  be  used  twice.  Inputs  are  R  and 
LAX  IS.  Returns  the  corresponding  rotation  sequence  ATT(1)  (0  to  2PI),  ATT(2)  (-PI 
to  PI),  ATT(3)  (0  to  2PI)  in  radians. 

RTOQ  This  subroutine  converts  a  nine  element  rotation  matrix  into  a  four  element  quaternion. 

RTPOSE  This  subroutine  transposes  the  rotation  matrix  R1  into  matrix  R2. 

VTOR  This  subroutine  converts  a  three  element  vector  into  a  nine  element  rotation  matrix  for 

coordinate  system  transformations.  The  NAXIS  axis  is  subjected  to  a  single  rotation 
which  points  it  in  the  direction  of  the  vector.  NAXIS  values  of  1, 2,  and  3  correspond 
to  the  x,  y,  and  z  axes  respectively.  Hie  resulting  rotation  matrix  R  is  returned.  R  may 
then  be  used  to  transform  to  the  new  system:  V  -  (R)(  V). 


4.2  M50TOT.FOR 

ASTRO  Calculates  M50  to  true-of-date  rotation  matrix  by  constructing  the  appropriate  precession 

and  nutation  rotation  matrices,  then  performing  matrix  multiplication. 

GMPRD  Multiplies  matrices. 

M50TOT  Subroutine  to  build  a  rotation  matrix  for  M50  to  true-of-date  coordinate  conversion. 

Input  four  digit  year,  month,  day,  and  GMT  (sec).  Returns  the  3x3  rotation  matrix 
called  RESULT. 

MJDT  Calculates  a  modified  Julian  day  number. 

NUTATE  Builds  the  nutation  rotation  matrix. 

PRECSS  Builds  the  precession  rotation  matrix. 

RMAT  Calculates  trigonometry  values,  then  builds  rotation  matrix. 

RMFU  Builds  rotation  matrix. 

SETCON  Initializes  the  math  constants  used  by  the  M50  converter. 
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5.  EXAMPLES 


5.1  Example  1. 

Convert  a  yaw,  pitch,  roll  sequence  to  the  corresponding  pitch,  yaw,  roll  sequence.  Show  that  the  two 
sequences  are  equivalent. 

This  example  is  illustrated  in  the  program  ATT.  Rotation  angles  are  given  values  in  a  data  statement 
and  are  stored  as  ATT!.  IAXIS1  holds  the  rotation  sequence  associated  with  ATT1.  The  values  3,2,1 
specify  the  2,y,x  axes  respectively.  This  corresponds  to  a  yaw,  pitch,  roll  sequence.  Similarly,  IAXIS2 
holds  the  rotation  sequence  associated  with  ATT2.  Values  of  2,3,1  correspond  to  a  pitch,  yaw,  roll 
sequence. 

After  changing  the  ATT1  angles  from  degrees  to  radians  we  are  ready  to  begin  the  conversion.  The 
attitude  sequence  defined  by  ATT1  and  IAXIS1  is  converted  to  a  rotation  matrix  using  a  call  to  ATOR. 
The  resulting  matrix,  TEMP,  is  then  converted  back  to  the  rotation  sequence  specified  by  IAXIS2  using 
the  call  to  RTOA.  The  values  of  interest  are  returned  in  ATT2. 

To  show  that  these  two  attitude  sequences  are  equivalent,  we  compare  them  using  the  logical  function 
AEQA.  A  value  of  true  or  false  is  returned  and  stored  as  RESULT.  Although  the  use  of  AEQA  is 
somewhat  trivial  in  this  example  program,  knowing  whether  or  not  two  attitude  sequences  are  equivalent 
is  a  common  need  and  so  it  is  demonstrated  here. 

The  rotation  sequences  and  their  associated  values  are  displayed  after  changing  back  from  radians  to 
degrees.  The  user  should  note  that  angles  and  axes  are  always  given  in  the  order  in  which  they  occur. 
It  is  the  values  of  IAXIS  which  contain  the  information  concerning  the  rotation  sequence. 

Finally,  the  equality  of  the  two  sequences  is  displayed  using  RESULT.  A  sample  output  follows  the 
program  listing. 

STS  MSID’s  of  interest  are  listed  in  Table  1. 
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ooo  o  oooo  oo  oooooo  noooo 


PROGRAM  ATT 


THXS  PROGRAM  CONVERTS  A  TAW,  PITCH,  ROLL  SSQUENCB  TO  A  PITCH,  YAW,  ROLL 
SEQUENCE.  XT  THEN  SHOWS  THAT  THE  TWO  SEQUENCES  ARB  EQUIVALENT. 


IMPLICIT  REAL*8(A-H,0~Z) 

LOGIC AIi  RESULT,  ASQA 

DIMENSION  ATT1(3),  IAXIS1(3),  ATT2(3),  XAXXS2(3),  TEMP(3,3) 

CHARACTER  LEGEND (3) *8 

DATA  LEGEND  /  'ROLL  »  ',  ’PITCH  «  ',  'YAW  «  •  / 

IAXXSl  VALUES  SPECIFY  YAW,  PITCH,  ROLL  SEQUENCE  FOR  ATTI. 

IAXIS2  VALUES  SPECIFY  PITCH,  YAW,  ROLL  SEQUENCE  FOR  ATT2. 

ATTI  DATA  IS  SUPPLIED  IN  DATA  STATEMENTS,  BUT  COULD  EASILY  BE  READ  FROM  A 
DATA  BASE. 

DATA  XAXIS1  /  3,  2,  1  / 

DATA  ATTI  /  0.3S82767D3,  0.2380823D0,  0.8965007D2  / 

DATA  XAXXS2  /  2,  3,  1  / 

CONVERT  DEGREES  TO  RADIANS  FIRST. 

TWOPX  -  8. DO  •  DATAN ( 1 . DO ) 

DTR  «  TWOPX  /  360. DO 
RTD  -  1.D0  /  DTR 
DO  100  1*1,3 

ATTI ( I )  «  ATTI ( X )  *  DTR 
100  CONTINUE 

CONVERT  ATTITUDE  SEQUENCE  1  TO  A  ROTATION  MATRIX,  THEN  CONVERT  THIS  TO 
ATTITUDE  SEQUENCE  2. 

CALL  ATOR(ATTl,IAXISl,TEMP) 

CALL  RTOA ( TEMP , ATT2 , XAXIS2 ) 

RESULT  -  AEQA ( ATTI , I AXIS1 , ATT2 , XAXXS2 ) 

DO  110  1-1,3 

ATTI ( I )  -  ATT1(I)  *  RTD 
ATT2  ( I )  -  ATT2  ( I )  *  RTD 
110  CONTINUE 

DISPLAY  RESULTS  TO  USER. 

WRITE (*,200)  (LEGEND ( IAXIS1 ( I ) ) , ATTI ( I) , X-l , 3) , 

+  (LEGEND ( XAXXS2 (I )),ATT2( I), 1-1,3) 

IF (RESULT)  THEN 

WRITE ( * , * )  '  EQUAL' 

ELSE 

WRITE ( * , * )  '  NOT  EQUAL' 

END  IF 

200  FORMAT (/,/, 2 ( IX, 3 ( 5X, A, F12 . 6) , / , / ) , / > 

END 


i 


YAW  -  358.276700 
PITCH  *  .238190 


PITCH  -  .238082  ROLL  -  89.650070 

YAW  -  >1.723285  ROLL  «  89.657233 


EQUAL 
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5.2  Example  2. 


Convert  M50  position  and  velocity  in  feet  per  second  to  Aries  true-of-date  position  and  velocity  in 
kilometers  per  second. 

This  example  is  illustrated  in  the  program  CNVM50.  In  this  example,  the  user  is  first  asked  to  specify 
input  and  output  file  names.  The  time,  position,  and  velocity  are  then  read  from  the  specified  input  file. 
A  printout  of  the  file  selected  for  this  example,  named  CNVM50.DAT,  is  included  at  the  end  of  the 
program  listing. 

The  time  of  day  is  converted  to  seconds,  then  a  cal!  is  made  to  M50TOT.  This  subroutine  uses  the  time 
information  to  calculate  rotation  matrix  values.  These  values  are  returned  as  the  3x3  array  M50TOD. 
The  user  may  be  interested  to  know  that  IMON,  IDAY  may  be  passed  as  either  1,  day  of  year  or  month, 
day  of  month, 

We  use  this  rotation  matrix  to  multiply  the  vector  of  interest  using  RMULTV.  In  this  way  the  M5G 
position  is  converted  to  TOD  position,  then  the  M50  velocity  to  TOD  velocity.  These  values  are 
converted  to  kilometers  and  kilometers  per  second  respectively,  and  displayed.  A  sample  display  is 
shown  following  the  program  and  input  file  listings. 

STS  MSID’s  of  interest  are  listed  in  Table  1. 


* 


* 
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oofs  non  oooo  ooo  ooo  ooo  oooo 


PROGRAM  CNVMSO 

THIS  PROGRAM  CONVERTS  POSITION  AND  VELOCITY  VICTORS  PROM  PT/8IC  IN  MSO  TO 
KM/SBC  IN  ARIES  TRUE  OP  DATI  COORDINATES. 


IMPLICIT  REAL *8  (A-H,0-Z) 

REALMS  M50TOD 

DIMENSION  POSM50 ( 3 ) ,  VELM50(3),  POSTOD(3),  VBLTOD ( 3 ) ,  M50TOD(3,3) 
CHARACTER* 12  IPILE,  OPILE 

FT/KM  CONVERSION 

DATA  PTKM  /  3280. 83 IDO  / 

USER  SELECTS  INPUT  AND  OUTPUT  PILE  NAMES 

WRITE (*,*)  '  ENTER  INPUT  PILB  NAME  • 

RBAD(*,  '  (A)  '  )  IPILE 

OPEN (UNIT- 1 , FILE-IFILB, STATUS- ' OLD ' , FORM* ' FORMATTED ' ) 

WRITE (*,*)  *  ENTER  OUTPUT  PILE  NAME  * 

READ(*, ’ (A) ’ )  OPILE 

OPEN ( UNIT-3 , FILE-OFILE , STATUS- ' NEW  * , FORM- • FORMATTED ’ ) 

READ  TIME  AND  VECTOR  INPUTS 

READ (1«*)  IYR, IMON, I DAY, IHR, IMIN/SBC 
READ ( 1 ,  * )  (POSMSO(I), 1-1,3) 

READ (lr*)  (VELMSO(I), 1-1,3) 

CONVERT  TIME  TO  SECONDS,  THEN  CONVERT  VECTORS  FROM  M50  TO  TRUE  USING 
ROTATION  MATRIX. 

GMT-IHR-3600 . DO+IMIN*  60 . DO+SEC 
CALL  M50TOT( IYR, IMON, IDAY,GMT,M50TOD) 

CALL  RMULTV(POSM50,M50TOD, POSTOD) 

CALL  RMULTV ( VELM50 , M50T0D , VBLTOD ) 

CONVERT  TO  XM 

DO  100  1-1,3 

POSTOD ( I )  -  POSTOD ( I )  /  PTKM 
VELTOD ( I )  -  VELTOD ( I )  /  PTKM 
100  CONTINUE 

WRITE  RESULTS 

WRITE (*,*)  POSTOD, VELTOD 

WRITE ( 3 , * )  IYR, IMON, IDAY , I HR, IMIN, SEC 

WRITE (3, * )  (POSTOD(I), 1-1,3) 

WRITE ( 3,  * )  (VELTOD(I), 1-1,3) 

CLOSE (UNIT-1) 

CLOSE (UNIT-3) 

END 


1990  1  306  8  IS 

0 . OOOOOOOOOOOOOOOB+OOO 

750. 624415039622700  5649.784942912215000  3398.785550422316000 

-5.275955101502905  -2.398301548788066  5.146722574414398 
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5.3  Example  3. 

Use  an  M50  to  body  axis  quaternion  and  body  axis  to  UVW  Euler  angles,  construct  the  M50  to  UVW 
rotation  matrix. 

This  example  is  illustrated  in  the  program  QM50.  The  quaternion  values,  Euler  angles,  and  axis 
sequence  information  were  taken  from  a  NASA  PATH  tape.  They  appear  as  data  statements  at  the 
beginning  of  the  program. 

The  M50  to  Body  quaternion  is  converted  to  a  rotation  matrix  using  QTOR.  Values  are  returned  as  the 
3x3  array  M50B.  After  changing  the  Euler  angles  from  degrees  to  radians,  the  Body  to  UVW  attitude 
sequence  is  converted  to  a  rotation  matrix  using  ATOR.  The  result  is  returned  as  the  3x3  array  BUVW. 
These  two  rotation  matrices  are  multiplied  using  RMULTR,  which  returns  the  3x3  array  M50UVW. 
This  last  rotation  matrix  is  the  result  of  two  successive  rotations,  M50B  followed  by  BUVW,  and  is  the 
desired  result. 

These  values  are  then  displayed  to  the  user.  A  sample  display  follows  the  program  listing.  The  reader 
might  be  interested  in  comparing  the  values  of  the  M50UVW  array  to  those  calculated  using  a  different 
method  and  shown  in  the  next  example. 

STS  MSID’s  of  interest  are  listed  in  Table  1. 
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non  o  no  o  no  oq 


PROGRAM  QM50 

C - - - 

C  THIS  PROGRAM  USES  A  QUATERNION  TO  CONSTRUCT  AN  M50  TO  BOOT  AXIS  ROTATION 
MATRIX.  IT  THEN  USES  EULER  ANGLES  TO  CONSTRUCT  A  BOOT  AXIS  TO  UVW  ROTATION 
MATRIX.  THESE  MATRICES  ARB  MULTIPLIED  TO  GXVE  AN  M50  TO  UVW  ROTATION 
MATRIX. 


IMPLICIT  REALMS (A-H,0-S) 

REAL *8  M50UVW,  M50B 

DIMENSION  QM50B<4),  M50B(3,3),  BUVW(3,3),  M50UVW(3,3),  IAXIS(3), 

+  ATT(3) 

FOR  THIS  EXAMPLE,  DATA  IS  SUPPLIED  IN  DATA  STATEMENTS.  IT  COULD  EASILY  BE 
READ  IN  FROM  A  DATA  BASE  INSTEAD.  SOURCE  OF  DATA  IS  NASA  PATH  TAPE 
51F-MISSION  TDRS  41D  WEST  -  SEG.  NO.  14.  IAXIS  VALUES  SPECIFY  YAW,  PITCH, 
ROLL  SEQUENCE  FOR  ATT. 

DATA  QM50B  /  0.2599793D0,  0. 54275S2D-1,  0.3427433D0,  >0.9011060  / 

DATA  IAXIS  /  3,  2,  1  / 

DATA  ATT  /  0.3582767D3,  0.2380823D0,  0.8965007D2  / 

CONVERT  DEGREES  TO  RADIANS  FIRST. 

TWOPI  ■  8. DO  *  DATAN(l.DO) 

DTR  -  TWOPI  /  360. DO 
DO  100  1-1,3 

ATT(I)  -  ATT(I)  <*  DTR 
100  CONTINUE 

CALL  QTOR(QM50B,M50B) 

CALL  ATOR(ATT, IAXIS, BUVW) 

CALL  RMULTR(M50B,BUVW,M50UVW) 

DISPLAY  RESULTS  TO  USER. 

WRITE (*,200)  'M50UVW  -  ' , ( (M50UVW( I, J) , J-l, 3) , 1-1,3) 

200  FORMAT (/ , /, IX, A, 3(T1S, 3 (5X, F10.6) ,/),/,/) 

END 


M50UVW  - 

-.844416 

.526901 

.096629 

-.282325 

-.591032 

.755628 

.455252 

.610783 

.647834 

t 
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5.4  Example  4. 

Use  M50  position  and  velocity  vectors  to  calculate  the  M50  to  UVW  rotation  matrix.  Use  Aries  true -of* 
date  position  and  velocity  vectors  to  calculate  the  True  to  UVW  rotation  matrix.  Combine  these  to 
obtain  the  M50  to  True  rotation  matrix,  then  check  independently. 

This  example  is  illustrated  in  the  program  UVW.  Values  for  M50  and  True  positions  and  velocities  were 
taken  from  a  NASA  PATH  tape.  They  appear  as  data  statements  at  the  beginning  of  the  program. 

The  M50  to  UVW  rotation  matrix  is  constructed  from  the  M50  position  and  velocity  vectors  using 
VTOR.  In  the  first  cal!  to  VTOR,  the  M50  coordinate  system  is  rotated  to  point  the  1,  or  x,  axis  in  the 
direction  of  the  position  vector.  This  rotation  makes  the  x'  axis  collinear  with  the  U  axis.  The  rotation 
matrix  necessary  for  this  operation  is  returned  as  the  3x3  array  Rl.  The  M50  velocity  vector  is  then 
transformed  to  the  new  system  using  RMULTV,  which  returns  the  transformed  velocity  as  TEMP.  This 
transformed  velocity  is  now  projected  onto  the  yV  plane  by  setting  the  x'  component  to  zero.  Next  the 
2,  or  y',  axis  is  pointed  in  the  direction  of  our  transformed,  projected  velocity  vector  by  using  VTOR  a 
second  time.  The  rotation  matrix  for  this  transformation  is  returned  as  the  3x3  array  R2.  This  second 
rotation  can  be  described  as  a  rotation  about  the  x*  axis  which  results  in  a  z"  velocity  component  of  zero. 
These  two  rotations  establish  conditions  sufficient  to  define  the  UVW  system.  The  rotation  matrices  Rt 
and  R2  are  multiplied  using  RMUI.TR.  The  product,  returned  as  the  3x3  array  M50UVW,  is  the  M50 
to  UVW  rotation  matrix  of  interest. 

The  true-of-date  to  UVW  rotation  matrix  is  constructed  from  True  position  and  velocity  vectors  in  a 
similar  manner,  and  returned  as  the  3x3  array  TODUVW.  Since  we  are  interested  in  going  from  UVW 
to  True,  RTPOSE  is  used  to  transpose  the  TODUVW  array.  The  result  is  returned  as  the  3x3  array 
UVWTOD. 

To  complete  the  construction  of  the  M50  to  True  rotation  matrix,  the  rotation  matrices  M50UVW  and 
UVWTOD  are  multiplied  using  RMULTR.  The  product  is  returned  as  the  3x3  array  M50TOD. 

Next,  the  M50  to  True  rotation  matrix  is  constructed  independently  using  the  date  and  time  associated 
with  the  position  and  velocity  vectors.  These  values  have  been  placed  in  a  data  statement  at  the 
beginning  of  the  program.  This  information  is  converted  to  a  rotation  matrix  using  M50TOT,  and  the 
results  are  returned  in  the  3x3  array  TEST. 

Finally,  we  check  to  see  if  the  two  rotation  matrices  are  equal  using  the  logical  function  REQR.  A  value 
of  true  or  false  is  returned  and  stored  as  RESULT.  The  values  of  interest  are  then  displayed  to  the  user. 
A  sample  display  follows  the  program  listing. 

STS  MSID’s  of  interest  are  listed  in  Table  1. 
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PROGRAM  UVW 

C - — - 

C  PROGRAM  TO  CALCULATE  M50  TO  UVW  AND  ARZI8  TRUB-OF-DATB  TO  UVW  ROTATION 

C  MATRICES  FROM  POSITION  AND  VELOCITY  VECTORS.  THE  TWO  MATRICES  ARB  THEN 

C  MULTIPLIED  TO  GIVE  A  M50  TO  TRUE  CONVERSION  MATRIX.  THIS  IS  CHECKED 

C  AGAINST  THE  M50TOT  SUBROUTINE  VALUES. 

C - - - - - 

c 

IMPLICIT  REAL*8(A-H,0-Z) 

LOGICAL  REQR,  RESULT 
REAL*8  K50UVW,  M50TOD 

DIMENSION  POSM50(3),  VBLM50(3),  POSTOD ( 3 ) ,  VELTOD (3),  TEMP(3) 

DIMENSION  Rl(3,3),  R2(3,3),  M50UVW(3,3),  TODUVW(3,3), 

+  UVWTOD (3,3),  MSOTOD (3,3),  TBST(3,3) 

C 

C  FOR  THIS  EXAMPLE,  DATA  IS  SUPPLIED  IN  DATA  STATEMENTS.  IT  COULD  EASILY  BE 

C  READ  FROM  A  DATA  BASE  INSTEAD.  SOURCE  OF  DATA  IS  NASA  PATH  TAPE 

C  SIP-MISSION  TORS  41D  WEST  -  SEG.  NO.  14. 

C 

DATA  IYR,  IMON,  IDAY,  GMTSBC  /  1985,  8,  1,  1001.86957DO  / 

C 

DATA  POSM50  /  -0. 5652093D4,  0.3526012D4,  0.6467874D3  / 

DATA  VELM50  /  -0.21788S3D1,  -0.45S3907D1,  0.5834034D1  / 

DATA  POSTOD  /  -0. 5681994D4,  0.3481981D4,  0.5274173D3  / 

DATA  VELTOD  /  -0.2162737D1,  -0.4581256D1,  0.5826428D1  / 

C 

C  CONSTRUCT  THE  M50  TO  UVW  ROTATION  MATRIX. 

C 

C  FIND  THE  ROTATION  MATRIX  WHICH  POINTS  THE  X  AXIS  IN  THE  DIRECTION  OF  THE 

C  M50  RADIUS  VECTOR.  THIS  TRANSFORMATION  DEFINES  THE  U  AXIS. 

CALL  VTOR ( POSM50 , Rl , 1 ) 

C  USE  THIS  MATRIX  TO  TRANSFORM  THE  VELOCITY  VEL'  -  (Rl)  (VEL) . 

CALL  RMULTV ( VELM50 , Rl , TEMP ) 

C  NOW  PROJECT  VEL*  ONTO  V,W  PLANE  BY  SETTING  THE  U  COMPONENT  TO  ZERO.  THEN 

C  FIND  ROTATION  MATRIX  WHICH  POINTS  THE  V  AXIS  IN  THE  DIRECTION  OF  THIS 

C  VELOCITY  PROJECTION.  THIS  IS  EQUIVALENT  TO  A  ROTATION  ABOUT  THE  U  AXIS 

C  WHICH  RESULTS  IN  NO  W  COMPONENT  OF  VELOCITY. 

TEMP ( 1 ) >0 • DO 

CALL  VTOR ( TEMP, R2, 2) 

C  COMBINE  THESE  TWO  ROTATIONS. 

CALL  RMULTR ( Rl , R2 , M50UVW ) 

C 

C  CONSTRUCT  THE  TOD  TO  UVW  ROTATION  MATRIX  IN  A  SIMILAR  MANNER. 

C 

CALL  VTOR ( POSTOD, Rl,l) 

CALL  RMULTV ( VELTOD, Rl, TEMP) 

TEMP ( 1 ) *0 • DO 

CALL  VTOR ( TEMP, R2, 2) 

CALL  RMULTR (Rl,R2,TODUVW) 

C 

C  NOW  CONSTRUCT  THE  M50  TO  TOD  ROTATION  MATRIX. 

C 

CALL  RTPOSE ( TODUVW , UVWTOD ) 

CALL  RMULrR(M50UVW, UVWTOD, M50TOD) 

C 

C  FIND  THIS  VALUE  INDEPENDENTLY  USING  SUBROUTINE  M50TOT. 

C 

CALL  M50TOT( IYR, IMON, IDAY, OMTSEC, TEST) 

C 

C  TEST  FOR  EQUALITY.  ACCEPTABLE  ERRORS  OF  ID-6  ARE  SPECIFIED  IN  LOOICAL 
C  FUNCTION  REQR. 

C 
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c 

c 

c 


RESULT  -  REQR ( TEST, M50T00 ) 
DISPLAY  RESULTS  TO  USER. 


( (M50UVW(I, J) , J-1,3) , 1-1,3) , 
( ( UVWTOD ( I, J) , J-l, 3) ,  1-1,3) , 
((M50TOD(I,J),J-1,3),I-1,3), 
<(  TBST( If J) , J-l, 3) , 1-1,3) 


WRITE (*,100)  'M50UVW  - 

+  ' UVWTOD  - 

+  'M50TOD  - 

+  'M50TOT  « 

IF (RESULT)  THEN 

WRITE ( * , * )  '  MSOTOD  EQUALS  M50TOT* 

ELSE 

WRITE;*,*)  •  MSOTOD  DOBS  NOT  EQUAL  M50TOT1 
END  IF 

100  FORMAT (/ , /,4 ( IX, A, 3 (T15, 3 ( 5X, F10. 6),/),/)) 
END 


M50UVW  - 

-.844416 

-.282325 

.455252 

.526901 

-.591032 

.610783 

.096629 

.755628 

.647834 

UVWTOD  - 

-.848883 

.520204 

.093735 

-.280239 

-.593280 

.754643 

.448179 

.614335 

.649406 

M50TOD  - 

.999963 

.007907 

.003438 

-.007907 

.999969 

.000018 

-.003437 

-.000045 

.999994 

M50TOT  - 

.999963 

.007907 

.003438 

-.007907 

.999969 

.000018 

-.003437 

-.000045 

.999994 

M50TOD  EQUALS  M50TOT 


t 


4 
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Table  1.  STS  Measurement/Stlmuli  Identification  for  Atlantis 


ID#  Description  Units 

V90W2310C  M50  TO  BODY  QUAT  TIME  S 

V90U2240C  M50-TO-MEASURED  BODY  QUAT  ELE  1  ND 

V90U2241C  M50-TO-MEASURED  BODY  QUAT  ELE  2  ND 

V90U2242C  M50-TO-MEASURED  BODY  QUAT  ELE  3  ND 

V90U2243C  M50-TO-MEASURED  BODY  QUAT  ELE  4  ND 

V90U2641C  M50  WRT  LVLH  QUAT  1  ND 

V90U2642C  M50  WRT  LVLH  QUAT  2  ND 

V90U2643C  M50  WRT  LVLH  QUAT  3  ND 

V90U2644C  M50  WRT  LVLH  QUAT  4  ND 

V97U2218C  ADI  REF  QUAT  ELEM  1  ND 

V97U2219C  ADI  REF  QUAT  ELEM  2  ND 

V97U2220C  ADI  REF  QUAT  ELEM  3  ND 

V97U2221C  ADI  REF  QUAT  ELEM  4  ND 

V98H1248C  BODY  WRT  INERTIAL  Q  WD  1  ND 

V98H1249C  BODY  WRT  INERTIAL  Q  WD  2  ND 

V98H1250C  BODY  WRT  INERTIAL  Q  WD  3  ND 

V98H1251C  BODY  WRT  INERTIAL  Q  WD  4  ND 

V90H2141C  BODY  ATTITUDE  ERROR-PITCH  DEO 

V90H2142C  BODY  ATTITUDE  ERROR-YAW  DEG 

V90H2143C  BODY  ATTHUDE  ERROR-ROLL  DEG 

V90H2202C  BODY  ROLL  ATTITUDE  EULER  ANGLE  DEG 

V90H2217C  BODY  PITCH  ATTITUDE  EULER  ANGLE  DEG 

V90H2230C  BODY  YAW  ATTITUDE  EULER  ANGLE  DEG 

V92H3333C  POR  P  ATT  ORB  STR  DEG 

V92H3334C  POR  YAW  ATT  ORB  STR  DEG 

V92H3335C  POR  R  ATT  ORB  STR  DEG 

V95H7467C  REQD  INERTIAL  ROLL  ANGLE  RAD 

V95H7468C  REQD  INERTIAL  PITCH  ANGLE  RAD 

V95H7473C  CURRENT  INERTIAL  ROLL  ANGLE  RAD 

V95H7474C  CURRENT  INERTIAL  PITCH  ANGLE  RAD 

V95H7475C  CURRENT  INERTIAL  YAW  ANGLE  RAD 

V95H7484C  REQD  INERTIAL  YAW  ANGLE  RAD 

V72R0916C  LH  ADI  ROLL  RATE  DEG/S 

V72R0917C  LH  ADI  PITCH  RATE  DEG/S 

V72R0918C  LH  ADI  YAW  RATE  DEG/S 

V72R1116C  AFT  ADI  ROLL  RATE  DEO/S 


Tablet,  (cont’d) 


ID#  Description  Units 

V72R1117C  AFT  ADI  PITCH  RATE  DEG/S 

V72R1118C  AFT  ADI  YAW  RATE  DEG/S 

V90R2223C  IMU  DERIVED  BODY  RATE  X-AXIS  DEG/S 

V90R2224C  IMU  DERIVED  BODY  RATE  Y-AXIS  DEG/S 

V90R2225C  IMU  DERIVED  BODY  RATE  Z-AXIS  DEG/S 

V92R3323C  ACT  POR  P  ROT  RATE  DEG/S 

V92R3324C  ACT  POR  YAW  ROT  RATE  DEG/S 

V92R3325C  ACT  POR  R  ROT  RATE  DEG/S 

V95R7476C  IMU  BODY  RATE  AROUND  X-AXIS  DEG/S 

V95R7477C  IMU  BODY  RATE  AROUND  Y-AXIS  DEG/S 

V95R7487C  IMU  BODY  RATE  AROUND  Z-AXIS  DEG/S 

V90H4277C  X-COMP  OF  FLTRS  CURR  POS  VCTR-TLM  FT 

V90H4278C  Y-COMP  OF  FLTRS  CURR  POS  VCTR-TLM  FT 

V90H4279C  Z-COMP  OF  FLTRS  CURR  POS  VCTR-TLM  FT 

V92H3417C  POR  X-POSITION  DISPLAY  IN 

V92H3418C  POR  Y-POSITION  DISPLAY  IN 

V92H3419C  POR  Z-POSITION  DISPLAY  IN 


V95H0185C  X-COMP  OF  CURRENT  SHUTTLE  POS  VCTR  FT 
V95H0186C  Y-COMP  OF  CURRENT  SHUTTLE  POS  VCTR  FT 
V95H0187C  Z-COMP  OF  CURRENT  SHUTTLE  POS  VCTR  FT 
V95H0862C  X-COMP  OF  CURRENT  TARGET  POS  VCTR  FT 
V95H0863C  Y-COMP  OF  CURRENT  TARGET  POS  VCTR  FT 
V95H0864C  Z-COMP  OF  CURRENT  TARGET  POS  VCTR  FT 


V90L2557C  SEL  TOTAL  X  VEL  M50  FT/S 

V90L2558C  SEL  TOTAL  Y  VEL  M50  FT/S 

V90L2559C  SEL  TOTAL  Z  VEL  M50  FT/S 

V92R3320C  ACT  POR  X  TRANS  RATE  FT/S 

V92R3321C  ACT  POR  Y  TRANS  RATE  FT/S 

V92R3322C  ACT  POR  Z  TRANS  RATE  FT/S 


V95L0190C  X-COMP  OF  CURRENT  SHUTTLE  VEL  VCTR  FT/S 
V95L0191C  Y-COMP  OF  CURRENT  SHUTTLE  VEL  VCTR  FT/S 
V95L0192C  Z-COMP  OF  CURRENT  SHUTTLE  VEL  VCTR  FT/S 
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